POLYHARMONIC  SMOOTHING  SPLINES  FOR  MULTI  DIMENSIONAL  SIGNALS  WITH 

l/|Mr-LIKE  SPECTRA 


Shai  Tirosh,  Dimitri  Van  De  Ville  and  Michael  Unser 

Biomedical  Imaging  Group 
Swiss  Federal  Institute  of  Technology  (EPFL) 

CH-1015  Lausanne,  Switzerland 

e-mail:  {Shai  .  Tirosh,  Dimitri  .VanDeVille,  Michael .  Unser}@epfl .  ch 


ABSTRACT 

Motivated  by  the  fractal-like  behavior  of  natural  images,  we 
propose  a  new  smoothing  technique  that  uses  a  regulariza¬ 
tion  functional  which  is  a  fractional  iterate  of  the  Laplacian. 
This  type  of  functional  has  previously  been  introduced  by 
Duchon  in  the  context  of  radial  basis  functions  (RBFs)  for 
the  approximation  of  non-uniform  data.  Here,  we  introduce 
a  new  solution  to  Duchon’s  smoothing  problem  in  multiple 
dimensions  using  non-separable  fractional  polyharmonic  B- 
splines.  The  smoothing  is  performed  in  the  Fourier  domain 
by  filtering,  thereby  making  the  algorithm  fast  enough  for 
most  multi-dimensional  real-time  applications. 

1.  INTRODUCTION 

In  many  signal  processing  systems,  denoising  algorithms 
are  applied.  Denoising,  or  smoothing,  which  are  qualita¬ 
tively  the  same  thing,  can  be  done  in  many  ways.  Tra¬ 
ditionally,  there  have  been  two  communities  dealing  with 
this  problem:  the  signal  processing  community  and  the 
statistics  community.  The  first  one  deploys  popular  algo¬ 
rithms  such  as  filtering  (e.g.,  Wiener  filtering)  and  wavelet 
thresholding  with  its  variations  [1],  The  statistics  commu¬ 
nity,  which  frequently  uses  non-uniform  grids,  adheres  to 
Bayesian  restoration  fe.g.  [2])  and  smoothing  splines  [3] 
(particularly  thin-plate  splines  [4]). 

In  fact,  one  can  think  about  denoising  as  taking  some  in¬ 
formation  out  of  a  signal.  The  information  we  want  to  take 
out  depends  on  a-priori  information  or  assumptions  about 
the  noiseless  signal.  A  common  assumption  is  that  the  en¬ 
ergy  proportion  of  the  signal  is  most  important  in  the  low 
frequencies,  while  the  high-frequency  components  contain 
mainly  noise. 

Multi-dimensional  data,  such  as  images,  usually  have 
high  local  correlation  in  all  directions.  Finding  a  suitable 
smoothing  algorithm  that  uses  this  fact  is  challenging.  The 
simplest  method  is  to  apply  a  separable  algorithm  (e.g.  by 
using  tensor  product  functions),  which  does  not  address  this 


correlation.  However,  their  implementation  tends  to  be  fast, 
and  the  mathematics  are  straightforward. 

Another  idea  is  to  use  RBFs,  such  as  Duchon’s  (to,  s) 
splines  [5],  providing  us  with  a  span  of  isotropic  power 
functions.  The  problem  is  formulated  in  variational  terms 
with  a  Tikhonov-like  regularizer.  However,  the  RBFs  are 
poorly  conditioned,  thus  making  it  very  difficult  to  imple¬ 
ment  when  there  are  many  data  points  as  is  typical  in  signal 
processing.  In  addition,  the  method  is  computationally  very 
expensive. 

Under  the  Gaussian  assumption,  there  is  a  well- 
known  equivalence  between  the  Bayesian  formulation  and 
regularization  techniques.  The  appropriate  regulariza¬ 
tion  for  fractal-like  signals — signals  with  a  spectra  like 
0(1/  |Mf) — is  a  fractional  iterate  of  the  Laplacian  which 
whitens  the  signal.  In  fact,  it  has  recently  been  demonstrated 
that  many  natural  images  have  this  kind  of  spectral  behav¬ 
ior  [6,7], 

Here,  we  propose  a  new  algorithm  to  solve  Duchon’s 
problem  on  a  uniform  grid  using  non-separable  fractional 
polyharmonic  B-splines  basis  functions.  These  functions 
are  localized  versions  of  RBFs,  and  thus  span  the  same 
space.  We  work  out  the  solution  for  the  fractional  case,  ob¬ 
taining  a  suitable  algorithm  for  fractal-like  signals. 

The  fact  that  we  are  dealing  with  a  uniform  grid  allows 
us  to  develop  a  fast  Fourier-based  filtering  algorithm.  In 
addition,  our  algorithm  can  work  in  any  number  of  dimen¬ 
sions  without  any  special  modifications.  Thus,  we  are  able 
to  solve  the  statistics  formulation  (Duchon’s  splines)  with 
signal-processing  techniques  (Fourier  filtering),  achieving  a 
fast  algorithm. 

2.  MATHEMATICAL  PRELIMINARIES 
2.1.  Polyharmonic  B-Splines 

Our  multi-dimensional  basis  functions  are  the  polyharmonic 
B-splines  [8-10].  We  denote  this  function  of  order  a  as 
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(a)  Order  1.5  (b)  Order  2.6 

Fig.  1.  Example  of  Polyharmonic  B-splines  in  2D.  Note  that  they  are  not  compactly  supported  and  decay  like 
0(\\x\\-d~2). 


Pa  ( x ),  and  specify  its  Fourier  transform: 


T  {Pa  (at)}  =  Pa  M 


||sin(^/2)|r 
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with  ll>  =  ,Ud)  G  where  sin(u;/2)  = 

(sin  (oq/2) ,  •  •  •  ,  sin  (uid/ 2)).  Note  that  the  order  of  those 
functions  can  also  be  fractional  (see  also  [11]).  Figure  1 
shows  an  example  of  two  polyharmonic  B-splines  in  2D. 

An  important  property  of  these  B-splines  is  the  convo¬ 
lution  relation  Pai  *  Pa2  =  Pa!+a2 .  which  follows  directly 
from  their  definition  in  the  Fourier  domain. 


2.2.  Fractional  Differentiation 

In  multiple  dimensions,  it  is  convenient  to  consider  the 
isotropic  fractional  differential  operator: 

dlf{x)< — >  IM|S/M  ,  (2) 


the  smoothing  operator,  /,  is  defined  by:  /  = 

argmin/6i2  {e^}  where: 

e2s=  (9  (Xi)  -  f  fa))2  +  A  ||<9*/  (x)\\l2  .  (5) 

a^jeS 

The  first,  signal  term  quantifies  the  distance  between  our  so¬ 
lution  and  the  given  measurements  g  (tc  j),  on  S  —  a  discrete 
localization  of  the  continuous  space,  i.e.  the  discrete  points 
at  which  we  are  observing  the  signal.  The  second,  smooth¬ 
ness  term  penalizes  the  lack  of  smoothness  of  the  solution. 
A  is  a  regularization  parameter,  making  a  balance  between 
the  two  terms  —  higher  A  means  more  smoothing  and  less 
fidelity  to  the  data,  and  vice  versa. 

Duchon  has  shown  that  the  solution  of  (5)  in  the  general 
case  is: 


/(*)  =  cikP  (x  —  Xi)  +  polynomial,  (6) 

XieS 


which  is  defined  in  the  sense  of  distributions.  The  discrete 
counterpart  is  the  (Laplacian-like)  finite  difference  operator: 


s/2 


A*  < — >  1 1 2  sin  (o->/2) || s  =  2s  (  ^sin(w,;/2)" 


(3) 
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It  is  easy  to  prove  the  following  formula  that  allows  us  to 
easily  differentiate  polyharmonic  B-splines: 


dlPa  ( X )  =  AlPa-s  (x)  .  (4) 


3.  SMOOTHING  FORMULATION 

Duchon’s  smoothing  formulation  is  variational  with  a 
Tikhonov-like  regularization  [5].  The  solution  of 


where  p[x)  is  the  Green  function  of  d%s:  d^sp(x)  =  8{x). 
In  the  Fourier  domain,  this  leads  to:  p[x)  < — >  /3(w)  = 
.  One  problem  with  p(ui)  is  its  singularity  at  zero  fre¬ 
quency.  In  the  frequency  domain,  this  function  looks  very 
similar  to  the  polyharmonic  B-splines  (except  for  the  sin 
term).  By  adding  the  sin  term,  the  function  is  tempered  at 
zero. 

We  can  show  that  these  functions  span  the  same  space 
iff  Xi  €  Zd,  which  is  mostly  the  case  in  signal  pro¬ 
cessing  applications  (It  is  easy  to  see  that  —  = 

Y^ke zd  d  (k)  P  (x  —  k)).  This  assumption  also  allows  us 
to  use  fast  Fourier  techniques.  Furthermore,  the  polyhar¬ 
monic  B-splines  of  order  2s  reproduce  polynomials  of  de¬ 
gree  [2s]  —1,  which  allows  us  to  also  take  care  of  the  second 
polynomial  term  in  (6). 


We  can  now  solve  Eq.  (5)  using  a  spline  representation 
for  our  signals: 

/(«)=  c(fe)$2  a(x~k).  (7) 

fee  zrf 

It  can  be  proven  that  this  choice  of  order  of  spline  gives 
a  globally  optimal  solution  among  all  possible  function 
spaces.  The  proof  is  a  generalization  of  the  one  in  [11], 

By  substituting  (7)  in  (5)  and  using  the  differentiation 
property  of  the  splines  (4)  we  obtain: 


Experimental  results  show  that  most  natural  images  fit  well 
to  this  model,  with  r  being  usually  between  1  to  1.5. 

Assuming  that  our  prior  signal  model  is  fractal,  we  can 
then  apply  a  fractional  derivative  of  power  s  =  t  to  whiten 
this  signal.  This  suggests  that  choosing  s  according  to  the 
fractal  power  of  a  signal  should  yield  better  results. 

This  value  is  calculated  by  taking  the  radial  frequency 
response  of  the  image,  transforming  it  into  log-log  coordi¬ 
nates  and  fitting  a  straight  line,  whose  slope  yields  r.  Since 
experimental  results  show  that  this  value  is  usually  frac¬ 
tional,  it  is  a  good  inducement  for  using  fractional  splines. 


e«  =  (9, 9 )  -  2  (g,  02s  *  c)  +  ( 02s  *  c,  02a  *  c) 

+A(AJ  *  c,  AJ  *  c  * /32s)-  (8) 


5.  RESULTS 


Taking  the  partial  derivative  with  respect  to  c  (partial  deriva¬ 
tive  of  a  vector)  and  equating  it  to  zero,  we  get: 

{(02s)'  *  9 )  (*)  =  ((02 s)'  *  (02 s)  *  c )  (at) 

+  (A  (A®)'  *  A*  *  (02s)'  *  c)  (at)  (9) 

Going  to  the  Fourier  domain,  we  find  the  solution: 

B2s  (e*w) 


/  M  = 


B2s  (&tw)  +  A  1 1 2 sin  (u>/2)\Y 


9  M,  (10) 


where  the  polyharmonic  B-spline  periodization  sequence 
Ij2s  ( ejU> )  is  the  Fourier  transform  of  a  sampled  polyhar- 
monic  spline.  It  can  be  expressed  as  the  periodized  version 
of  (32s  (u>)  (using  Poisson  summation  formula): 


B2s  (e>u) 


=  fi2s  + 2nk ) 

Ate  zd 


=  E 

fee  zd 


1 1  sin  (tu/2  +  7rfc)||2s 
||cu/2  +  7rfc|j2s 
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We  demonstrate  our  approach  on  a  slice  of  an  MRI T2*  vol¬ 
ume,  as  used  in  functional  analysis  of  brain  activity  (this 
image  is  in  fact  the  brain  of  one  of  the  former  Ph.D  students 
of  our  group).  The  original  image  is  given  in  Fig.  2(a).  In 
Fig.  2(c),  we  show  the  results  of  smoothing  the  noisy  image 
of  Fig.  2(b).  In  this  case,  s  was  estimated  from  the  spectrum 
of  the  image,  and  we  chose  the  value  of  A  that  achieves  the 
best  SNR. 

There  are  known  methods  for  selecting  the  value  of  A, 
but  it  is  beyond  the  scope  of  this  paper.  Note  that,  when  A  is 
small,  the  result  fits  the  observation  more  closely,  and  edges 
are  preserved  better;  however,  this  comes  at  the  expense  of 
more  residual  noise.  On  the  other  hand,  when  A  is  too  large, 
the  image  is  over  smoothed,  and  edges  are  not  well  pre¬ 
served;  however,  we  suppress  most  of  the  noise.  The  right 
balance  will  depend  on  the  particular  application. 

Figure  3(a)  shows  the  radial  frequency  response  of  the 
image  in  Fig.  2,  and  the  regressed  fractal  model  (on  a  log- 
log  scale).  Note  that  the  value  of  s  for  the  smoothing  was 
chosen  according  to  this  analysis.  In  Fig.  3(b),  we  show  that 
this  way  of  selecting  s  is  indeed  optimal. 


This  sequence  can  be  also  regarded  as  the  autocorrelation 
sequence  of  0S.  We  have  found  an  efficient  way  to  calculate 
it,  by  applying  a  two-scale  relation  in  the  Fourier  domain. 
(More  details  about  this  will  be  published  in  a  future  paper). 

4.  FRACTAL-LIKE  BEHAVIOR  OF  NATURAL 
IMAGES 

Natural  phenomena  are  widely  presumed  to  be  self-similar 
from  a  statistical  perspective.  As  a  consequence,  natural  im¬ 
ages  tend  to  be  scale  invariant  —  seeing  an  object  from  two 
yards  or  one  yard  will  result  in  very  similar  images  trans¬ 
mitted  by  our  visual  system  [7], 

The  notion  of  self-similarity  is  reminiscent  of  fractals, 
which  possess  this  property.  This  leads  to  a  fractal-like 
model  for  natural  images  —  a  spectral  density  function  that 
behaves  like  O  (l/||u;||T),  where  r  is  the  fractal  order  [6], 


6.  CONCLUSIONS 

We  proposed  a  multi-dimensional  smoothing  algorithm  us¬ 
ing  non-separable  fractional  polyharmonic  B-splines.  This 
method  is  a  solution  of  Duchon’s  smoothing  formula¬ 
tion,  and  combines  the  statistics  community’s  methods 
with  the  signal  processing  community’s  tools,  allowing  for 
fast  smoothing  of  multi-dimensional  signals  using  Fourier- 
based  filtering.  This  is  only  possible  because  we  are  work¬ 
ing  on  a  uniform  grid  and  because  we  found  a  fast  algorithm 
to  compute  (11). 

Secondly,  we  use  the  fact  that  natural  images  behave  in  a 
fractal-like  way  with  respect  to  their  radial  spectrum,  which 
in  turn  allows  us  to  tune  the  value  of  s  —  the  fractional 
derivative  order  in  the  variational  smoothing  formulation. 
Additionally,  preliminary  experimental  results  confirm  that 
this  choice  is  optimal.  Our  approach  is  especially  suitable 


(a)  Original  image  (b)  Noisy  image  (c)  Optimal  least  squares  smoothing 

Fig.  2.  Smoothing  results.  SNR:  (b)  14.7804  (c)  22.4707,  with  A  =  3.07,  s  =  1.46633. 


Approx,  radial  freq.  behavior  of  image  and  regressed  line 


Distance  from  origin,  a  =  -1 .46633  ;  b  =  0.614891  (log  scale) 


(a)  Radial  frequency  behavior  (log-log  scale) 


Best  SNR  per  s  value 


(b)  Best  SNR  vs.  value  of  s 


Fig.  3.  Fractal-like  properties,  (a)  The  radial  frequency  behavior  of  the  image,  with  the  corresponding  regressed  line  (on  a 
log-log  scale),  (b)  The  effect  of  choosing  the  correct  value  of  s  on  the  optimal  denoising  SNR. 


for  fractal-like  signals  as  well  as  real  world  images  thanks 
to  our  fractional  extension  of  smoothing  splines. 
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